Cosmogenic nuclide dating of two stacked ice masses: Ong Valley, Antarctica

We collected a debris-rich ice core from a buried ice mass in Ong Valley, located in the Transantarctic Mountains in Antarctica. We measured cosmogenic nuclide concentrations in quartz obtained from the ice core to determine the age of the buried ice mass and infer the processes responsible for the emplacement of the debris currently overlaying the ice. Such ice masses are valuable archives of paleoclimate proxies; however, the preservation of ice beyond 800 kyr is rare, and therefore much effort has been recently focused on finding ice that is older than 1 Myr. In Ong Valley, the large, buried ice mass has been previously dated at > 1.1 Ma. Here we provide a forward model that predicts the accumulation of the cosmic-ray-produced nuclides ¹⁰Be, ²¹Ne, and ²⁶Al in quartz in the englacial and supraglacial debris and compare the model predictions to measured nuclide concentrations in order to further constrain the age. Large downcore variation in measured cosmogenic nuclide concentrations suggests that the englacial debris is sourced both from subglacially derived material and recycled paleo-surface debris that has experienced surface exposure prior to entrainment. We find that the upper section of the ice core is 2.95 + 0.18 $/$ −0.22 Myr old. The average ice sublimation rate during this time period is 22.86 + 0.10 $/$ −0.09 m Myr⁻¹, and the surface erosion rate of the debris is 0.206 + 0.013 $/$ −0.017 m Myr⁻¹. Burial dating of the recycled paleo-surface debris suggests that the lower section of the ice core belongs to a separate, older ice mass which we estimate to be 4.3–5.1 Myr old. The ages of these two stacked, separate ice masses can be directly related to glacial advances of the Antarctic ice sheet and potentially coincide with two major global glaciations during the early and late Pliocene epoch when global temperatures and CO₂ were higher than present. These ancient ice masses represent new opportunities for gathering ancient climate information.

1 Introduction

Ice cores from glaciers and ice sheets are used as an archive for paleoclimate proxies, including atmospheric gases, chemical compounds, and airborne particles (Dansgaard et al., 1969; Fredskild and Wagner, 1974; Castellano et al., 2004; Willerslev et al., 2007); however, the potential age of ice core records is limited by the fact that ice sheets are subject to deformation, ice flow, and basal melting. The oldest ice that has been recovered from the thickest parts of the Antarctic ice sheets is 800 000 years old (Jouzel et al., 2007). Although it is hypothesized that ice up to 1–1.5 Ma may also exist at great depth in the ice sheet (Fischer et al., 2013), recovering this ice would be a complex and costly endeavor. Therefore, we currently lack archives of climate information that extend beyond ∼ 0.8 Ma.

Bare ice is, in general, thermodynamically unstable under typical atmospheric pressure–temperature conditions and therefore prone to melt and/or sublimate. However, there exist regions of topographically constrained, extremely slow ice flow in which ice up to 2.7 Ma has been recovered near the surface (Yan et al., 2019). There are also several areas within the Transantarctic Mountains where glacial ice is covered by supraglacial debris. A thick debris cover thermally insulates the ice surface and provides a physically plausible means of preserving near-surface ice for long periods. For example, Sugden et al. (1995) found glacier ice in Beacon Valley underlying a supraglacial debris containing 8.1 Myr old volcanic ash, leading them to conclude that the ice is older than 8.1 Ma. However, the antiquity of the Beacon Valley ice has been questioned on the basis of data suggesting that ice lost to sublimation is dynamically replaced by ice flow from upstream glaciers, resulting in a situation where relatively young ice underlies relatively old debris (Van Der Wateren and Hindmarsh, 1995; Ng et al., 2005; Hindmarsh et al., 1998; Stone et al., 2000). The lack of ice older than ∼ 1 Ma severely limits our direct paleoclimate record and creates uncertainties when modeling future climate predictions which include modeled configuration of the past Antarctica Ice Sheet (Bulthuis et al., 2019; Noble et al., 2020). This is particularly important during the Pliocene epoch (Dolan et al., 2018; Haywood et al., 2009), in which global surface temperatures and CO₂ levels were higher than present (Pagani et al., 2010; Seki et al., 2010) and which is considered analog for current anthropogenic warming.

In Ong Valley, Miller Range, Transantarctic Mountains, a mass of glacier ice at least several tens of meters thick is found buried underneath < 1 m of supraglacial debris. Cosmogenic nuclide measurements from the supraglacial debris suggest an age of > 1.1 Ma but most likely > 1.8 Ma (Bibby et al., 2016). We collected a 944 cm long ice core from this buried ice mass and use concentrations of ¹⁰Be, ²¹Ne, and ²⁶Al from the englacial debris to further constrain the age of the ice mass.

Our goal is to determine the age of the ice, understand its overall geologic history, and evaluate its potential use as a paleoclimate archive. We present a novel dating application of cosmogenic nuclides which aims to quantify a complex exposure history of this buried ice mass. By comparing measured cosmogenic nuclide concentrations from the englacial and supraglacial debris with modeled concentrations, the nuclide inventory inherited from prior exposure can be distinguished from that produced after ice emplacement. We then apply a cosmogenic nuclide burial dating method to the inherited inventory as an age constraint. We show that two sections of the ice core contain recycled surface debris that can be burial dated. The upper section is ∼ 3 Ma, which we interpret as the emplacement age of the bulk of the buried ice. The lower section has a significantly older burial age of > 4 Ma, and we interpret it as a portion of an older ice mass either in situ or transported during emplacement of the younger ice.

2 Study area

Ong Valley is a ∼ 1.5 km wide and ∼ 7 km long glacial valley located in the Miller Range of the central Transantarctic Mountains, Antarctica (83.25^∘ S, 157.72^∘ E). The current valley floor gradually rises from an elevation of 1500 m above sea level (m a.s.l.) to 1700 m a.s.l. at the valley head. Over the span of 1 year (2011), the recorded air temperature in the valley ranged between −49.0 and −4.0 ^∘C, with a mean of −23.9 ^∘C (Bibby et al., 2016). In the head of the valley is a small alpine glacier, and the valley mouth is blocked by a 2 km wide exposed glacial ice front of the Argosy Glacier (Fig. 1).

Figure 1 Location of Ong Valley. (a) Cropped USGS 1:250 000 scale topographic map of Miller Range, Antarctica, showing the location of Ong Valley. The red rectangle indicates the location of Ong Valley opening perpendicular to Argosy Glacier. (b) WorldView 2 satellite image of Ong Valley, Antarctica (© 2016 Maxar). The dots indicate sampling sites for pit and ice core (orange), Middle drift surface boulders (cyan), and lateral moraine boulders (blue). The legend shows the prefixes of sample names.

The valley floor is mostly covered by a well-developed and distinctive system of three glacial drifts, referred to as Young, Middle, and Old (Bibby et al., 2016) (Figs. 1 and 2). These deposits were first described in 1975 and later identified as soil chronosequences, increasing in age and maturity with distance from the Argosy Glacier (Mayewski, 1975; Scarrow et al., 2014). Bibby et al. (2016) found that the three drift units were ablation tills formed by sublimation of debris-rich glacier ice that advanced into the valley. Eventually the ice became stagnant and began to sublimate, which led the englacial debris to accumulate on the surface as supraglacial debris. Although some of the supraglacial debris in Ong Valley could originate from a rockfall or colluvium from adjacent slopes, the drifts either have convex topography (Middle and Young drifts) or are bounded by prominent moraine ridges (Old drift), and therefore significant input from local slopes is only possible immediately adjacent to valley walls. In addition, surfaces of active glaciers in the region uniformly lack significant surface sediment. While eolian sediment transport onto the drifts is possible, drift surfaces are mainly composed of clasts and boulders too large for eolian transport.

Figure 2 Oblique aerial photograph of Ong Valley (left) with added markings (right) indicating the Young, Middle, and Old drifts (Bibby et al., 2016). The photograph is looking northward and down-valley. The dots represent the sampling sites as identified in Fig. 1.

Exposure dating of the supraglacial debris from each drift has revealed the ages of 11–13 ka (Young), > 1.1 Ma (Middle), and > 2.7 Ma (Old) (Bibby et al., 2016). The Young and Middle drifts have buried ice under 0.1–0.5 m (Young) and 0.6–0.8 m (Middle) of loose supraglacial debris completely concealing the ice. In this paper we refer to the buried ice below the Middle drift as the Middle ice. In contrast, the Old drift is devoid of buried ice, which, presumably, has sublimated over extended exposure (Bibby et al., 2016). The highest surface elevation of these buried ice masses in the valley is located on top of the Middle drift and ∼ 200 m above the current Argosy Glacier surface elevation at the valley mouth. Bibby et al. (2016) used cosmogenic nuclide data from the surface debris layer to estimate that sublimation rates of 19–23 m Myr⁻¹ and surface erosion rates of 0.7–0.9 m Myr⁻¹ have persisted in the valley since deposition of the drifts. All three drifts have related lateral moraines on the valley walls that trace the original elevation of the ice surface. The Old drift also has a distinct end moraine close to the head of the valley, which shows no signs of influence from the small alpine glacier currently located at the head of the valley.

Major and minor mineral analyses of englacial debris and the supraglacial debris from each of the three drifts display shared provenance that are significantly different from the local bedrock of Ong valley (Edwards et al., 2014; Morgan et al., 2020). The local bedrock is primarily dominated by Hope Granite and the Argosy Gneiss (Barrett et al., 1970). This indicates that these drifts were deposited by past advances of the Argosy Glacier.

3 Cosmogenic nuclide applications relevant for dating Ong Valley buried ice

We use measured concentrations of cosmic-ray-produced nuclides both in the supraglacial debris covering the buried ice mass in Ong Valley and in debris within the buried ice in order to determine the age of the ice, its sublimation and erosion rates, and the geologic history of englacial debris. Cosmogenic nuclides are rare isotopes produced by cosmic-ray interactions with matter at Earth's surface (Dunai, 2010). Because the cosmic-ray flux is rapidly attenuated with depth below the surface, the concentration of cosmogenic nuclides can be used to date or quantify geologic processes that form or bury surface materials (Lal, 1991; Nishiizumi et al., 1993). For example, nuclide concentrations in rocks or sediment brought rapidly from depth to the surface and not subsequently disturbed can be interpreted as an “exposure age” for rock or sediment. If such material is then re-buried, the decay of cosmogenic radionuclides can be used to compute a “burial age”. In general, given a set of assumptions derived from the geologic context of a sample, measured cosmogenic nuclide concentrations can yield information about the timing of geologic events or the rates of geologic processes such as erosion and sedimentation that have affected that sample.

There are three means by which we can apply cosmogenic nuclide data to determine the age of the buried ice in Ong Valley. First, we can apply exposure dating to the supraglacial debris; this approach was taken by Bibby et al. (2016). However, because the supraglacial debris is formed by sublimation after the deposition of glacial ice, the duration of the surficial exposure of the supraglacial debris is, by definition, shorter than the age of the ice. Thus, in this case, a surficial exposure dating is expected to yield only a minimum age. We can exposure date the supraglacial debris by determining the apparent exposure age of the surface, defined as the exposure age calculated from a nuclide concentration with the assumption of a single period of exposure continuing until today without erosion or burial. The basic assumptions for such exposure dating are that a sample (i) has never been exposed to cosmic rays prior to entrainment; (ii) has only been exposed to cosmic rays since deposition; and (iii) has then never been covered, displaced, nor disturbed while exposed at the surface. Further, we can use two nuclides to quantify both erosion and exposure time (Lal, 1991). The concentration of each nuclide is a function of exposure age and surface erosion rate. Therefore, measurements of both nuclides yield two independent equations in which both unknowns can be solved for in certain circumstances.

The second approach that we use to determine the age of the buried ice mass is generally referred to as “depth-profile dating” (Hidy et al., 2010), and it involves measurements of both surface and subsurface nuclide concentrations. This approach relies on the observation that surface concentrations show a greater dependence on sublimation and/or surface erosion rates and lesser dependence on emplacement age compared to concentrations in the subsurface below several meters depth (Stone et al., 1998; Braucher et al., 2009). Thus, paired subsurface and surface measurements can, in principle, yield a unique solution for both ice age and erosion rate. In Ong Valley, we apply this approach to both the supraglacial debris and subsurface englacial debris in core by creating a forward model that predicts nuclide concentrations at all depths as a function of the age of a deposit, the surface erosion rate of the deposit, and, in this case, the ice sublimation rate leading to formation of the supraglacial debris deposit. Fitting this forward model to a data set then yields best-fitting estimates of these input parameters.

The third approach can be used if any of the englacial debris has formerly been exposed at the Earth's surface and subsequently buried. Then we can apply a burial dating method based on the decay of cosmogenic radionuclides produced during the initial period of exposure. The principle of burial dating is that different cosmogenic nuclides, such as ¹⁰Be and ²⁶Al, in the same mineral are produced at a fixed ratio that depends on the production rates. Samples that have experienced a single period of exposure at the surface have nuclide concentration ratios corresponding to the production ratio. If the two nuclides have different half-lives, burial of the sample to a depth at which the cosmic-ray flux is diminished causes the observed ratio to change through time due to the different rates of radioactive decay (Lal, 1991). Thus, given several other assumptions, the ratio reflects the duration of burial. Although we had no prior reason to expect the englacial debris to be sourced from recycled surface debris, we show later that it is, in fact, the case in Ong Valley. Therefore, at our field site, burial dating can be used as an approach to constrain the age of the ice.

The general organization of this paper follows the process of starting with the simplest expectation that is supported by the geology and learning and adapting as the results are unfolding in the course of the project. Therefore, we start with the simplest assumption that the deposit of the Middle ice is from a single glacial advance. We first measure the cosmogenic nuclide concentrations in debris obtained from the ice core. We then provide a forward model that predicts the accumulation of cosmogenic nuclides at depth below the surface. This forward model includes the processes responsible for the emplacement of the supraglacial debris. Fitting of the forward model to the cosmogenic nuclide data set reveals that the deepest section of measured nuclide concentrations is incompatible with an assumption of a single glacial advance. This leads us to apply additional analysis that was not originally planned for but which is compatible with the two mapped glacial advances into the upper valley. Therefore, the application of burial dating of an older ice layer is introduced in a latter section as an addition to the forward model. This addition allows us to model and explain the observed downcore variations in the cosmogenic nuclide concentrations which are supported by the local geology.

4 Methods

4.1 Sample collection

During the Austral summer, 2017/2018, we collected (i) pit samples from unconsolidated supraglacial debris at drill site 17-OD1, (ii) an ice core taken directly below the pit samples, and (iii) erratic boulders from other locations on the Middle drift surface and correlative lateral moraines (Fig. 1). The drill site OD1 was located at a central high point within the Middle drift. We chose this site because any deformation of the buried ice should be minimized at this location, and colluvium and rockfall from the valley walls cannot reach the site. We determined the location and elevation of the core site using postprocessed differential GPS. Boulder samples were located using uncorrected handheld GPS, and their elevations were checked against the Reference Elevation Model of Antarctica (REMA) digital elevation model (DEM) (Howat et al., 2019). Topographic shielding calculations for the sites follow Balco et al. (2008, with accompanying online material).

4.1.1 Vertical pit sampling of surficial regolith

We excavated a hand dug pit for sampling the vertical section of the supraglacial debris (Fig. 3a). The pit was located in the center of a patterned ground polygon formation in which the surface did not show signs of reworking caused by former active polygon boundaries. The supraglacial debris is a sandy diamict with clasts of all sizes up to large boulders. The debris surface is covered by a lag deposit of clasts larger than approximately 5–10 cm. The clasts are mostly angular with occasional faceted and/or weakly polished surfaces. Clast lithologies include both local bedrock and other rock types not locally present. The supraglacial debris shows no sign of stratification nor presence of ice-cemented regolith, and a sharp boundary can be observed between the debris and the underlying debris-rich ice mass (Fig. 3a).

Figure 3 Photographs of the drill site 17-OD1. (a) Excavated pit for sampling of the supraglacial debris. The Middle drift ice surface is found at the bottom of the pit. The yellow ruler in the image measure ∼ 60 cm from the bottom of the pit to the surface of the supraglacial debris. (b) Winkie drill installed above the excavated pit for ice coring. Photograph looking south.

Bulk sediment samples spanning 3–8 cm in sample depth were collected at approximately 10 cm depth intervals (Table 1). A total of six samples were collected throughout the pit, from the surface to the ice–debris boundary. The depth of the pit is 62 cm reaching from the surface of the supraglacial debris down to the ice surface.

4.1.2 Ice core sampling (and visual observations)

We collected a 944 cm long and 7 cm diameter vertical ice core directly below the deepest supraglacial debris sample. This allows us to construct a depth profile of cosmogenic nuclide concentrations extending from the surface of the supraglacial debris to the bottom of the ice core. An updated Winkie drill (Fig. 3b) was used for drilling in mixed-media ice (debris-rich ice) (Boeckmann et al., 2020). The recovered ice core was divided into segment lengths of 21 cm or less and stored frozen in individual watertight sample containers at the drill site. The samples were kept frozen at <− 20 ^∘C until processed in the geochemistry laboratory at the University of North Dakota.

4.1.3 Sampling of Middle drift boulders

Samples of glacially transported boulders were collected during the 2011/2012 field season (Middle drift boulders) and the 2017/2018 field season (lateral moraine boulders) (Fig. 4). Boulders were selected for sampling based on their evidence for stability in the drift. We preferentially selected boulders that were partially buried in supraglacial debris and that otherwise showed no signs of overturning due to cryoturbation processes post-deposition. We also selected boulders that rose > 50 cm from the debris surface to limit burial by snow. For the Middle drift samples (sample prefixes 11-OV-ER^*, Fig. 1), we identified samples in the center of the valley to avoid any rockfall from the valley wall. For the lateral moraine boulders (sample prefixes 17-OV-ERR^*, Fig. 1) on the valley wall, we avoided rectangular boulders that were perched on the debris surface as these appeared to be jointed boulders more recently eroded out of the bedrock. At all sites, deeply weathered boulders were avoided in an attempt to limit samples with complex exposure histories. All boulders sampled were either granite or gneiss that contained a high percentage of quartz. With a hammer and chisel, we removed a 1–2 kg sample from the top of each boulder. GPS location, exposure geometry, and sample thickness were recorded. Boulders in the middle of the drift should provide minimum exposure ages for the Middle drift as they were exposed as the Middle ice sublimated.

4.2 Ice core visual observations

Visual inspection of the ice core indicates that it is primarily composed of debris-rich ice in which enclaves, bands, and pods of sediment-rich ice, as well as individual mineral grains, are clearly visible (Figs. 5 and S1 in the Supplement). The englacial debris is poorly sorted, ranging from clay size to clasts exceeding the diameter of the borehole (> 7 cm). The core variously includes isolated mineral grains within relatively clear ice, centimeter-scale enclaves of moderately sorted sandy sediment, and centimeter-scale enclaves of very poorly sorted clay-rich diamict resembling glacial debris as typically found in wet-based glacial settings.

Most sections of debris-rich ice are banded, with the banding arising from variations in sediment concentration and with a thickness of the individual bands ranging from a few millimeters to centimeters. The orientations of these debris bands are variable but commonly steeply dipping. The drill system did not maintain core orientation during core extraction, and therefore, the dip direction is unknown. Although the debris concentration and sorting were highly variable, we did not observe any systematic downcore variation in visually identifiable sediment properties such as color or mineral composition. The ice core also contains three 1.2–1.5 m segments of mainly clear ice between 130 and 250, 410 and 560, and 710 and 840 cm below the surface. Although these segments include a few isolated debris bands and pods, they are primarily composed of clear ice lacking visually identifiable englacial debris.

Many of these features, particularly the discontinuous, variable, and steeply dipping nature of the debris bands, are characteristic of ice near the bed of glaciers and ice sheets where subglacial- or englacial-derived debris is inhomogeneously mixed into the ice by shear deformation. Basal ice at the present margin of the Argosy Glacier in lower Ong Valley has similar characteristics. However, the currently exposed section of the ice front displays a lower density of debris bands and a higher density of clear ice than the ice core. Therefore, we interpret the majority of the ice core as being consistent with characteristics of basal ice expected during an advance of the Argosy Glacier into central Ong Valley.

Figure 4 Photographs of surface boulders sampled from (a–c) the Middle drift surface and (d–f) the east lateral moraine; (a) 11-OV-ER-117, (b) 11-OV-ER-118, (c) 11-OV-ER-119, (d) 17-OV-ERR-213, (e) 17-OV-ERR-217, and (f) 17-OV-ERR-218.

4.3 Sample processing

4.3.1 Sample preparation

Each ice core segment containing visible englacial debris was weighed before the ice was melted at room temperature. The liquid was then separated from the englacial debris, and the debris was oven dried at ∼ 70 ^∘C. Then the samples were re-weighed. Resulting masses of ice and debris in each core segment were then used to compute the density of the core using Eq. (4). The ice core alternates between sections of clean ice and debris-rich ice with a maximum debris concentration of 0.57 by weight (Fig. 5). Ice core segments with little or no debris were assumed to have the density of ice, 0.917 g cm⁻³.

Figure 5 Image of the 17-OD1 ice core. The full image is stitched together from multiple individual pictures each covering approximately 20 cm of core length. The blue graph shows the corresponding debris concentration calculated by weight at depth below the surface. The topmost 62 cm is the thickness of the overlying supraglacial debris layer (not shown in pictures). See the complete ice core image in the Supplement.

The density of the supraglacial debris was measured by packing the sediment sample into a measuring cup of known volume and weighing it. This process was repeated five times for each sample to obtain a representative mean density of the supraglacial debris of ∼ 1.8 g cm⁻³, excluding grains > 2 cm.

For cosmogenic nuclide analyses, the pit samples from the supraglacial debris were sieved to separate a grain size fraction of 250–500 µm, and the englacial debris was sieved to a grain size of 250–833 µm. This was done to maximize the amount of sample for cosmogenic nuclide analyses. Quartz grains were isolated and cleaned following the procedure described in Stone (2004). When needed, adjacent segments were merged so that enough quartz was available to permit precise Be and Al measurements (see Table 1 for details).

4.3.2 Cosmogenic nuclide extraction and analysis

After quartz preparation but prior to dissolution for Be–Al extraction, a ∼ 0.5 g aliquot of the prepared quartz was split off for ²¹Ne measurements. These employed the “Ohio” noble gas mass spectrometer and extraction line at the Berkeley Geochronology Center (BGC). Details of neon (Ne) isotope measurements on this system are described in Balco and Shuster (2009) and Balter-Kennedy et al. (2020). ²¹Ne concentrations in replicate analyses of the CRONUS-A intercomparison standard (Vermeesch et al., 2015) measured during analytical sessions in this study ranged from 314.3 ± 9.4 to 320.8 ± 6.1 × 10⁶ atoms g⁻¹, indistinguishable from the accepted value of 320 ×10⁶ atoms g⁻¹. In Tables 1 and 2 and in the Supplement, we report ²¹Ne concentrations as excess ²¹Ne relative to atmospheric composition. Excess ²¹Ne includes both cosmogenic ²¹Ne and, potentially, nucleogenic ²¹Ne derived from uranium (U) and thorium (Th) decay. U and Th concentrations in quartz from Ong Valley lithologies (Sams, 2016), expected Ne closure ages of these lithologies, and the observation that ²¹Ne concentrations in drift and boulder samples from the Young drift are significantly higher than expected for the last glacial maximum (LGM) age of that drift all indicate that nucleogenic ²¹Ne is significant in quartz in these lithologies. For samples from supraglacial and englacial debris used to fit our forward model for nuclide accumulation, we did not make a correction for nucleogenic ²¹Ne because it would be equivalent to inherited cosmogenic ²¹Ne in our model simulations. Thus, such a correction would not affect values for ages or process rates inferred from the model simulations. In calculating exposure ages and erosion rates for boulder samples, we corrected for nucleogenic ²¹Ne using an estimate of 7 ± 3 ×10⁶ atoms g⁻¹ obtained from ²¹Ne data on boulders of similar lithology from the Young drift (Sams, 2016). With this estimate, nucleogenic ²¹Ne comprises 2 %–7 % of total excess ²¹Ne in boulder samples. In addition, as discussed below, we applied this correction to subsurface samples used for burial dating.

Chemical extraction and preparation of beryllium and aluminum from remaining quartz extracted from drill core samples were performed at the University of Vermont/National Science Foundation Community Cosmogenic Facility following the process described in Corbett et al. (2016). The pit and ice core samples were processed in two separate batches of 12 samples, in which each batch included a process blank and a standard. For all samples, 250 µg ⁹Be was added with a beryl carrier made at the facility with a concentration of 291 µg mL⁻¹. In addition, a ²⁷Al carrier commercially available as an inductively coupled plasma (ICP) standard from SPEX with a concentration of 1000 µg mL⁻¹ was added only to samples having < 1500 µg of total Al. The amount of ²⁷Al carrier added was based on the total amount of native ²⁷Al in a sample quantified by inductively coupled plasma optical emission spectroscopy (ICP-OES) analysis. Quartz isolation and beryllium extraction for boulder samples followed the same procedure, except that three of the boulder samples (11-OV-ER-117, 118, 119) were processed in chemistry laboratories at the Center for Accelerator Mass Spectrometry, Lawrence Livermore National Laboratory (LLNL-CAMS).

Ratios of ¹⁰Be $/$ ⁹Be measured at Lawrence Livermore National Laboratory (LLNL) are normalized to the 07KNSTD3110 standard (Nishiizumi et al., 2007) with an assumed ¹⁰Be $/$ ⁹Be ratio of 2.85 × 10⁻¹². All ²⁶Al $/$ ²⁷Al ratios were measured at PRIME and normalized to the KNSTD-01-5-2 standard (Nishiizumi, 2004) with an assumed ²⁶Al $/$ ²⁷Al ratio of 1.818 × 10⁻¹².

Both ¹⁰Be and ²⁶Al measurements were corrected for background using a procedural blank measured in each batch. Procedural blanks were run with samples from the core site and were 14.9 ± 3.3 × 10³ and 49.8 ± 9.4 × 10³ atoms ¹⁰Be and 259 ± 45 × 10³ and 37 ± 20 × 10³ atoms ²⁶Al, and blanks run with boulder samples ranged from 24 to 109 × 10³ atoms ¹⁰Be. For the majority of samples, blank corrections account for less than 1 % of total ¹⁰Be or ²⁶Al atoms present. The exception is several samples of englacial debris from the ice core, for which blank corrections were up to 4 % of total ²⁶Al atoms present and up to 9 % of total ¹⁰Be atoms present. The reported uncertainty in the measured nuclide concentrations accounts for all sources of analytical errors, including AMS measurement uncertainties, concentration measurement of ¹⁰Be and ²⁶Al, and procedural blanks. Measurement details appear in the Supplement.

4.4 Forward exposure model

The exposure history of the Middle drift in Ong Valley is complex, and the nuclide production cannot be accounted for by simply exposure dating the supraglacial debris. Therefore, we apply a forward model which attempts to account for the geological processes that result in the thickening of the supraglacial debris and accumulation of cosmogenic nuclides at depth. As described in Bibby et al. (2016), the concentration of cosmogenic nuclides in the ice mass and the supraglacial debris is expected to result from a series of events: (i) debris-rich glacial ice was deposited into Ong Valley during glacial advancement; (ii) the ice mass became stagnant and began to sublimate, which caused the englacial debris to accumulate on the ice surface as a supraglacial debris layer; (iii) as the ice continued to sublimate additional debris was added to the supraglacial debris layer from below, bringing deeper samples closer to the surface; (iv) at the same time the supraglacial debris layer was subjected to surface erosion at a rate slower than the accumulation of debris from sublimation such that the supraglacial debris thickness increased with time. The present-day thickness of the supraglacial debris layer is therefore a function of the age of ice emplacement, rate of ice sublimation, concentration of debris in ice, and rate of surface erosion. The numerical forward model attempts to account for the series of events listed above which lead to the cosmogenic nuclide concentrations measured today at depth below the surface in the Middle ice. Input parameters to the model include the age of ice emplacement, the sublimation rate of the ice, and the surface erosion rate of the supraglacial debris layer. The debris concentration in the ice is not an independent input parameter but is determined from the age of the ice, sublimation and erosion rates, and present-day thickness of the supraglacial debris (Sect. 4.4.2). The model then predicts nuclide concentrations in supraglacial and englacial debris. By fitting the model to the observed nuclide concentrations, we obtain estimates for the age of the ice and for sublimation and erosion rates.

4.4.1 Shielding mass

The accumulation of cosmogenic nuclides in a target mineral at and below the surface is dependent on the cosmic-ray flux. The cosmic-ray flux is significantly attenuated as it travels through mass, such as ice and mineral matter, which leads to the decreased production rate of cosmogenic nuclides below the surface. Therefore, the production rate of cosmogenic nuclides at and below the surface is dependent on the amount of shielding mass above the given sample. The shielding mass is dictated by the amount of debris in the ice, as well as the densities of the debris and ice, and, hence, becomes the foundation of the forward model. This section describes the calculation of the shielding mass. The resulting production rate is introduced in a later subsection as it is based on the shielding mass which is temporally variable.

The shielding mass is the cumulative mass of sediment and ice overlying each sample per unit area and has units of grams per square centimeter (g cm⁻²). It is equal to the product of the sample depth (cm) and the mean density of the overlying material (g cm⁻³). The density is related to the concentration of the suspended debris C_D and the ice C_I in each core segment and is calculated from the total segment weight, M_t (g), and dried sediment weight, M_s (g), such that

$$\begin{array}{}\text{(1)}& {\displaystyle }{C}_{\mathrm{D}}={\displaystyle \frac{M_{\mathrm{s}}}{M_{\mathrm{t}}}},\text{(2)}& {\displaystyle }{C}_{\mathrm{I}}=\mathrm{1}-C_{\mathrm{D}}.\end{array}$$

The density of the ice core can then be calculated by mixing the two ice core components based on volume. Assuming the density of the ice ρ_I to be 0.917 g cm⁻³ and the density of the debris ρ_D to be that of rock, 2.68 g cm⁻³, for an ice core segment weight of 1 g, the total weight of the ice core segment is then

$$\begin{array}{}\text{(3)}& M_{\mathrm{t}}=C_{\mathrm{D}}+C_{\mathrm{I}},\end{array}$$

resulting in the density of the ice–debris mixture in an ice core segment ρ_M being

$$\begin{array}{}\text{(4)}& {\mathit{\rho }}_{\mathrm{M}}={\displaystyle \frac{\mathrm{1}}{\left(\frac{C_{\mathrm{D}}}{{\mathit{\rho }}_{\mathrm{D}}}\right)+\left(\frac{C_{\mathrm{I}}}{{\mathit{\rho }}_{\mathrm{I}}}\right)}}.\end{array}$$

4.4.2 Depth as a function of time

Once an ice mass with a mixture of ice and debris is emplaced at some time in the past defined as T (year), the shielding mass above a given sample in the ice is decreasing through time at the rate determined by the sum of the sublimation rate of the ice and the surface erosion rate of the supraglacial debris. As the sample reaches the top of the ice, it becomes part of the supraglacial debris and then approaches the surface solely at the rate of surface erosion. The sublimation rate, s (cm yr⁻¹), is defined as a constant rate in which the surface of the ice–debris mixture (bottom of the supraglacial debris layer) is lowering. Note that this parameter represents a surface lowering rate due to sublimation and is not the same as a sublimation rate of pure ice as would be considered in a thermodynamics context. Therefore, the initial surface of the ice mass is sT (cm) above the present surface.

The rate at which mass is being lost by sublimation is the product of the sublimation rate and the density of the sublimating material. Since the ice mass consists of a mixture of ice and debris, only part of the ice mass is sublimating. The rate of mass loss associated with sublimation is given by s(1−C_D)ρ_M (g cm⁻² yr⁻¹). While the ice is sublimating, the debris suspended in the ice mass is left behind on the ice surface and accumulates as supraglacial debris. The rate at which mass is added to the bottom of the supraglacial debris by sublimation is then sC_Dρ_M (g cm⁻² yr⁻¹). By assuming that at the time of emplacement, the thickness of the supraglacial debris above the ice mass was zero, then with constant sublimation rate and erosion rate the total mass thickness of the supraglacial debris, Z_till (g cm⁻²), created by ice sublimation can be expressed as

$$\begin{array}{}\text{(5)}& Z_{\mathrm{till}}=Ts{C}_{\mathrm{D}}{\mathit{\rho }}_{\mathrm{M}}-TE,\end{array}$$

where T is the age in years before present that the ice was emplaced, and E is the erosion rate expressed in mass units (g cm⁻² yr⁻¹). Equation (5) leads to the constraint that sC_Dρ_M>E as the thickness of the supraglacial debris cannot be negative.

From field measurements, the thickness of the supraglacial debris is known to be 110 g cm⁻² for 17-OD1-Pit2. Therefore, for any arbitrary values of age, sublimation rate, and erosion rate, the debris concentration must be chosen such that the measured supraglacial debris mass thickness is obtained after T years. Assuming the ice mass mixture only consists of ice and debris, then the term C_D and ρ_M are different representations of the debris mass embedded in the sublimating ice.

Multiplying Eq. (4) by C_D, the debris concentration of the lost mass associated with sublimation can be solved.

$$\begin{array}{}\text{(6)}& C_{\mathrm{D}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{D}}\left(C_{\mathrm{D}}{\mathit{\rho }}_{\mathrm{M}}\right)}{{\mathit{\rho }}_{\mathrm{D}}{\mathit{\rho }}_{\mathrm{I}}-\left({\mathit{\rho }}_{\mathrm{I}}-{\mathit{\rho }}_{\mathrm{D}}\right)\left(C_{\mathrm{D}}{\mathit{\rho }}_{\mathrm{M}}\right)}}\end{array}$$

Equations (5) and (6) are two independent equations including the term C_Dρ_M. By isolating and substituting the term C_Dρ_M in Eq. (5) into Eq. (6) the debris concentration now becomes independent of the density of the ice core and is instead a function of sublimation rate (s), erosion rate (E), and age of ice emplacement (T).

$$\begin{array}{}\text{(7)}& C_{\mathrm{D}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{D}}\left(\frac{Z_{\mathrm{till}}+TE}{Ts}\right)}{{\mathit{\rho }}_{\mathrm{D}}{\mathit{\rho }}_{\mathrm{I}}-\left({\mathit{\rho }}_{\mathrm{I}}-{\mathit{\rho }}_{\mathrm{D}}\right)\left(\frac{Z_{\mathrm{till}}+TE}{Ts}\right)}}\end{array}$$

The debris concentration is constrained such that $\mathrm{0}\ge C_{\mathrm{D}}\le \mathrm{1}$. Further, Eq. (5) leads to the constraint that sC_Dρ_M>E as the thickness of the supraglacial debris cannot be negative.

When predicting the concentration of cosmogenic nuclides, it is crucial to know a sample's depth at present, defined as Z_S,now in units of mass depth (g cm⁻²), and its depth at some time, t (yr), in the past, Z(t). Since the samples collected consist of both ice core samples and sediment from the above-lying supraglacial debris, there are two separate cases of how a sample has approached the surface in the past.

In case 1, the sample is in the ice at present such that the sample depth is greater than the depth of the supraglacial debris, Z_s,now > Z_till. From time of ice emplacement, the sample has then approached the surface at the rate of the ice sublimating and the rate of surface erosion such that

$$\begin{array}{}\text{(8)}& z\left(t\right)=Z_{\mathrm{s},\mathrm{now}}+ts\left(\mathrm{1}-C_{\mathrm{D}}\right){\mathit{\rho }}_{\mathrm{M}}+tE.\end{array}$$

In case 2, the sample is in the supraglacial debris at present such that Z_s,now < Z_till. The time that the sample has been in the supraglacial debris is then defined as the mass height of the sample above the supraglacial debris base depth divided by the rate of mass addition to the supraglacial debris.

$$\begin{array}{}\text{(9)}& T_{\mathrm{till}}={\displaystyle \frac{\left(Z_{\mathrm{till}}-Z_{\mathrm{s},\mathrm{now}}\right)}{s{C}_{\mathrm{D}}{\mathit{\rho }}_{\mathrm{M}}}}\end{array}$$

By the time of emplacement, the sample has then approached the surface in the same way as above. However, once the sample reaches the top of the ice and becomes part of the supraglacial debris, it approaches the surface at the rate of surface erosion only. This further allows for two scenarios to occur.

In case 2a, if at time t the sample is in the supraglacial debris such that t< T_till, then

$$\begin{array}{}\text{(10)}& z\left(t\right)=Z_{\mathrm{s},\mathrm{now}}+tE.\end{array}$$

In case 2b, if at time t, the sample is in the ice, where t> T_till, then

$$\begin{array}{}\text{(11)}& z\left(t\right)=Z_{\mathrm{s},\mathrm{now}}+T_{\mathrm{till}}E+\left(t-T_{\mathrm{till}}\right)\left(s\left(\mathrm{1}-C_{\mathrm{D}}\right){\mathit{\rho }}_{\mathrm{M}}+E\right).\end{array}$$

With this, it is now possible to calculate the depth of a sample at any given time in the past since the ice emplacement. This is needed in order to predict the total accumulation of nuclides in a sample, which is dependent on the nuclide production at depth.

4.4.3 Cosmogenic nuclide production at depth

The cosmogenic nuclides ¹⁰Be, ²⁶Al, and ²¹Ne are produced by high-energy spallation, negative muon capture, and fast muon interactions (Dunai, 2010). The high-energy spallation particles are likely to react with mass in the atmosphere and at Earth's surface. Therefore, the production of cosmogenic nuclides due to spallation reaction is highest at the surface and considerably decreases with depth. Muons are much less likely to interact with mass and therefore travel farther below the surface before stopping (Lal, 1991). While cosmogenic nuclide production at depth below the surface is solely due to muon production, muons are responsible for less than 1 % of the total production at the surface for all nuclides (Balco, 2017).

We calculated the ¹⁰Be production rate using the “LSDn” scaling method (Lifton et al., 2014) as implemented in version 3 of the online exposure age calculator originally described by Balco et al. (2008) and subsequently updated, as well as the CRONUS-Earth “primary” calibration data set (Borchers et al., 2016). This yields a long-term (> 1 Myr) average production rate due to spallation of 25.3 atoms ¹⁰Be g${}_{\mathrm{qtz}}^{-\mathrm{1}}$ yr⁻¹ (where g_qtz referes to gram of quartz). We then assumed that the ²¹Ne $/$ ¹⁰Be production ratio is 4.03 (Balco et al., 2019) and the ²⁶Al $/$ ¹⁰Be production ratio is 6.75 (Balco et al., 2008).

The spallation production rate at the surface P_sp(0) for a given cosmogenic nuclide decreases exponentially with depth (Lal, 1991) such that

$$\begin{array}{}\text{(12)}& P_{\mathrm{sp}}\left(z\right)=P_{\mathrm{sp}}\left(\mathrm{0}\right)e^{\left(\frac{-z}{\mathrm{\Lambda }}\right)},\end{array}$$

where z is the mass depth (g cm⁻²), and L is the attenuation length, defined as the distance in which the energetic cosmic-ray flux intensity reduces by a factor of e⁻¹ due to scattering and absorption processes. The attenuation value varies depending on altitude and latitude and is taken to be 140 g cm⁻² in Antarctica for ¹⁰Be, ²⁶Al, and ²¹Ne (Borchers et al., 2016; Balco et al., 2019).

Muon production rates are not well approximated by a single exponential function. As depth increases, the energy of the remaining muons that have not yet stopped is higher, and therefore it takes proportionally longer for those to stop. The calculations of the production rates due to negative muon capture follows that of Heisinger et al. (2002a) and production rates due to fast muon interactions according to Heisinger et al. (2002b), and they are combined into a total muon production rate at depth, P_μ(z). The surface topography surrounding the sampling site also shields the samples of cosmic rays and will need to be accounted for when computing the production rate. This topographic shielding scaling factor S_G is applied to the spallation and not the muon production rate (Balco et al., 2008). A topographic shielding of 0.993 was measured for drill site 17-OD1. The total production rate as function of depth can then be described as

$$\begin{array}{}\text{(13)}& P\left(z\right)=S_{\mathrm{G}}{P}_{\mathrm{sp}}\left(z\right)+P_{\mathit{\mu }}\left(z\right).\end{array}$$

In Fig. 6 we show the calculated changes in mass depth and production rate for a sample collected in the supraglacial debris 50 cm below the surface. The following arbitrary but illustrative model parameter values are used for ice emplacement age, sublimation rates, and erosion rates: 1 Ma, 20 m Myr⁻¹, and 0.1 m Myr⁻¹, respectively. The supraglacial debris thickness is that measured at drill site 17-OD1 and is 62 cm (110.15 g cm⁻²). From the time of ice emplacement, a sample's depth has decreased linearly due to ice sublimation and surface erosion as the age of the ice increases, with a distinct change in rate once the sample exits the ice and becomes part of the supraglacial debris following Eqs. (10) and (11). It is also observed that such samples experience great changes in nuclide production rates, following Eq. (13), from time of emplacement until present.

Figure 6 Graphical representation of the temporal change in mass depth (blue line) and production rates (orange lines) since ice emplacement for a sample found today at 90 g cm⁻² (∼ 50 cm below the surface) in the supraglacial debris. Here, the given sample has an initial shielding mass of 1600 g cm⁻² (∼ 14 m) and approaches the surface at the combined rate of sublimation and erosion (Eq. 11). Once the sample reaches the ice surface it becomes part of the supraglacial debris (∼ 0.2 Ma) and approaches the surface solely at the rate of surface erosion (Eq. 10). As the sample's mass depth decreases, it experiences considerable changes in production rates of cosmogenic nuclides (Eq. 13). The arbitrary model parameter values for these calculations are 1 Ma ice emplacement, 20 m Myr⁻¹ sublimation rate, and 0.1 m Myr⁻¹ erosion rate.

4.4.4 Predicting cosmogenic nuclide concentration at depth

When exposed to cosmic rays, a sample will begin to accumulate cosmogenic nuclides over exposure time T such that the total accumulation, N (atoms g⁻¹), in a subsurface sample at depth z (g cm⁻²) can be expressed as the integral of the production rate a sample undergoes from time of ice emplacement to the present. Since the exposure history of the englacial debris goes beyond the exposure history of the ice mass, some amount of inherited background nuclides N_inh (atoms g⁻¹) are present. While the concentration of the stable cosmogenic nuclide ²¹Ne continues to build up, some of the unstable radionuclides, ¹⁰Be and ²⁶Al, are lost to radioactive decay. This is expressed as an exponential such that the total number of cosmogenic radionuclides for a sample can be calculated using Eq. (14) and simplified to Eq. (15) for the stable nuclide ²¹Ne.

$$\begin{array}{}\text{(14)}& {\displaystyle }{N}_i=N_{i,\mathrm{inh}}{e}^{\left(-{\mathit{\lambda }}_{i}T\right)}+{\int }_{\mathrm{0}}^T{P}_i\left(z\left(t\right)\right)e^{\left(-{\mathit{\lambda }}_{i}t\right)}\mathrm{d}t\text{(15)}& {\displaystyle }{N}_{\mathrm{21}}=N_{\mathrm{21},\mathrm{inh}}+{\int }_{\mathrm{0}}^T{P}_{\mathrm{21}}\left(z\left(t\right)\right)\mathrm{d}t\end{array}$$

The subscript i refers to the radionuclide of interest, ¹⁰Be or ²⁶Al, and λ is the decay constant for the radionuclide i. The decay constants used in this paper are 4.99 × 10⁻⁷ and 9.83 ×  10⁻⁷ for ¹⁰Be and ²⁷Al, respectively. We evaluate these integrals numerically using the default algorithm (integral) in MATLAB.

Given a set of environmental conditions – (i) age of ice emplacement, (ii) sublimation rate of ice, (iii) surface erosion rate of the supraglacial debris, and (iv) inherited nuclide concentration for each of the nuclides – as well as measurable site conditions (e.g., elevation and supraglacial debris thickness), Eqs. (14) and (15) predict the current total cosmogenic nuclide concentration in a sample at a unique depth z at present time.

4.4.5 Model fitting statistics

We defined a misfit statistic to compare observed nuclide concentrations with those predicted by the model as the reduced χ² statistic weighted by the relative uncertainty in measurements of all three nuclide concentrations in each sample. A best fit is found using the constrained nonlinear multivariable optimizing function (fmincon) in MATLAB while optimizing for the free parameters: (i) the age that the ice was emplaced, (ii) sublimation rate of the ice since emplacement, (iii) surface erosion rate of the accumulating supraglacial debris, and (iv) the inherited nuclide concentrations in the englacial debris at the time of ice emplacement.

Uncertainty distributions on the best-fit values of the model parameters are derived from a 10 000-iteration Monte Carlo simulation. Each Monte Carlo iteration draws from a set of normally distributed uncertainty for the measured nuclide concentrations and uses an initial guess for the free model parameters based on published values for the Middle drift in Ong Valley (Bibby et al., 2016): 1.83 Ma ice emplacement age, 22.7 m Myr⁻¹ sublimation rate, 0.89 m Myr⁻¹ surface erosion rate, and inherited nuclide concentrations of 0.14 × 10⁶ atom g${}_{\mathrm{qtz}}^{-\mathrm{1}}$, 11.4 × 10⁶ atom g${}_{\mathrm{qtz}}^{-\mathrm{1}}$, and 0.82 × 10⁶ atom g${}_{\mathrm{qtz}}^{-\mathrm{1}}$ for ¹⁰Be, ²¹Ne, and ²⁶Al, respectively. Although uncertainties in calculated nuclide production rates are similar in magnitude to the uncertainty measured in nuclide concentrations, they are not included in these results as all the samples are from the same location and therefore must have the same surface production rate. Further, any of the Monte Carlo simulation results which neither converged nor satisfied the optimization function evaluation or resulted in an ice emplacement age younger than the last glacial maximum (LGM) (< 20 ka) were excluded.

5 Results

5.1 Measured ice core nuclide concentrations

We measured cosmogenic nuclide ¹⁰Be, ²¹Ne, and ²⁶Al in 6 supraglacial debris pit samples (17-OD1-PIT2^*) and 12 ice core samples (17-OD1-C1^*). In addition, ²¹Ne was measured for the surface sample (17-OD1-surf). The results show that cosmogenic nuclide concentrations in the englacial debris show large downcore variations (Table 1). The relative magnitude of the variation is larger for ¹⁰Be and ²¹Ne than for ²⁶Al, which indicates that the amount of time that has elapsed since ice emplacement is most likely on the order of several half-lives of ²⁶Al (the ²⁶Al half-life is 0.7 Myr).

The measured cosmogenic nuclide concentrations for ¹⁰Be, ²¹Ne, and ²⁶Al decrease monotonically from the surface samples to 412 cm below the surface, with the exception of three consecutive samples (17-OD1-C1-70-100, 17-OD1-C1-107-125, and 17-OD1-C1-125-145) that exhibit anomalously high nuclide concentrations between 132 and 207 cm below the surface (Table 1). At 562 cm below the surface the cosmogenic nuclide concentrations either remain constant or increase towards the bottom end of the core, dependent on the nuclide. The measured ²⁶Al concentrations of the englacial debris show some of the lowest concentrations measured from the ice core and remain constant at approximately 0.45 × 10⁶ atoms g${}_{\mathrm{qtz}}^{-\mathrm{1}}$. In contrast, the measured ¹⁰Be and ²¹Ne concentrations increase towards the bottom of the core and display much higher concentrations that are similar to those measured in the top 1 m of the ice core (¹⁰Be) and lower section of the pit (²¹Ne) in the supraglacial debris.

Table 1 Measured cosmogenic nuclide concentrations in quartz extracted from supraglacial debris (prefix 17-OD1-surf/PIT2^*) and ice core (prefix 17-OD1-C1^*).

5.2 Boulder surface exposure results

We measured cosmogenic nuclide ¹⁰Be and ²¹Ne in quartz from erratic boulders from the Middle drift surface and east lateral moraine (Fig. 1). Apparent ¹⁰Be and ²¹Ne exposure ages for the Middle drift surface boulder samples are 1.55–2.16 and 0.82–1.39 Myr for the east lateral moraine boulders (Table 2). Further, the apparent ¹⁰Be exposure ages reported for all boulder samples appear younger than of those from ²¹Ne and is an indication that some process (e.g., erosion, burial, etc.) must have occurred which decreases the ¹⁰Be nuclide concentration relative to ²¹Ne

Table 2 Exposure age data for boulders on the surface of Ong Valley Middle drift and correlative lateral moraines. DD signifies decimal degree.

Notes: (1) rock density is assumed to be 2.57 g cm⁻³ based on measurements on like lithologies in Ong Valley. (2) Exposure ages are calculated using default production rates and “LSDn” scaling in version 3 of the online exposure age calculator described by Balco et al. (2008) and subsequently updated. (3) Both internal (including only measurement uncertainty) and external (in parentheses; also includes production rate uncertainty) uncertainties are shown for apparent exposure ages. (4) ²¹Ne exposure age calculations include subtraction of 7 ± 3 ×10⁶ atoms g⁻¹ non-cosmogenic ²¹Ne (see text).

6 Model fitting

In the following sections we first highlight several important features of the data from the ice core and supraglacial debris that we seek to explain using the forward model for nuclide accumulation described above. We then fit the forward model to the observations and thereby obtain estimates for the emplacement age and sublimation rate of the buried ice. Lastly, we calculate a minimum age for the Middle drift in Ong Valley.

6.1 Qualitative observations of the ice core data

The forward model described in Sect. 4.4.4 is an exposure model that is based on an assumption that the englacial debris is well mixed and therefore all the samples contain the same amount of inherited nuclides. The model calculates the post-depositional nuclide production during the exposure (Eqs. 14 and 15). This is compatible with an ice mass that is emplaced during single glacial advance. After a single event of exposure, in which sublimation and erosion have occurred, the concentration for a given cosmogenic nuclide in the englacial debris must decrease monotonically with depth as the production rate decreases with increased shielding mass.

However, multiple sections of the ice core show an increase in cosmogenic nuclide concentration at depth (Table 1) and are therefore not compatible with a single exposure history described above. The most likely explanation for the observed increases in cosmogenic nuclide concentrations with depth in OD1 samples is that englacial debris in various sections of the core have variable exposure histories prior to entrainment in the ice and therefore have different inherited nuclide concentrations. Based on the measured nuclide concentrations we make the following two observations that guide our forward modeling.

The first observation is that the set of samples that display monotonically decreasing nuclide concentrations across the entire profile (Fig. 7, sections S1, E1, and E2) are segments of the ice core with relatively low debris concentrations. These samples follow the expectation of having an exposure history as outlined in the forward model. We define these samples to be of a low-nuclide concentration that is largely composed of subglacially derived material sourced from upstream of the Argosy Glacier and to have minimal surface exposure prior to entrainment. The debris from such samples is therefore identified to be of “englacial debris” (samples denoted englacial E1 and E2 and highlighted in blue in Fig. 7) including the current surface debris layer (S1) that we assume has also originated from englacial debris.

Figure 7 Depth plot of (a) sediment concentration (weight sediment / weight sediment and ice) and (b) measured ¹⁰Be, ²⁶Al, and ²¹Ne nuclide concentrations. The bold vertical lines represent the average measured cosmogenic nuclide concentrations in the samples. The width of the surrounding black boxes shows the measured uncertainty in the cosmogenic nuclide concentrations, and the horizontal black lines indicate the original ice core segments and are the same for all three nuclides. Color shades highlight the source of the debris. Yellow depicts the current supraglacial debris samples (S1) and has same origin as blue samples (E1 and E2) identified as “englacial debris” having no prior surface exposure before being subglacially entrained. Samples in shades of red are identified as “paleo-surfaces” (S2, S3, and S4) and are made of recycled surface debris that was exposed on the surface and entrained as the glacier overrode previously ice-free areas as it advanced into Ong Valley.

The second observation is that the debris-rich sections of the ice core (Fig. 7a, sections S2, S3, and S4) have much higher nuclide concentrations and do not conform to the assumptions of a well-mixed ice mass containing subglacially derived debris and only experiencing in situ englacial accumulation of cosmogenic nuclides. We consider the debris in these sections to be from one or more high-concentration sources that are likely composed of material that was exposed for an extended period of time on the surface and entrained as the glacier overrode previously ice-free areas as it advanced into Ong Valley. In fact, these samples have ²¹Ne concentrations similar to modern surface material (S1) and therefore must have been exposed for millions of years at or near the surface. Therefore, debris from these subsurface samples having higher nuclide concentrations is identified as “paleo-surfaces” and is observed in three distinct units. We use the term “paleo-surface” in reference to samples that are assumed to have previously been exposed at or near the surface prior to entrainment by the Middle ice. Our usage of the term does not imply that the surface is still in situ. Although these units must all contain some fraction of recycled surface debris, this may be of variable origin, so we identify these as S2, S3, and S4 (highlighted in shades of red in Fig. 7). In this paper we refer to this debris as recycled paleo-surface.

Based on the above observations we can only explain a subset of the data with the forward model. However, the presence of significant and variable inheritance in the debris from different sources suggests that we can apply a burial constraint to the inherited nuclide inventory having prior exposure. When fitting the model to the post-depositional nuclide inventory, this allows for an additional constraint on the age of ice emplacement.

6.2 Forward modeling used to explain the data set

The forward model predicts the accumulation of cosmogenic nuclides ¹⁰Be, ²⁶Al, and ²¹Ne in the ice core and the overlaying supraglacial debris during a single event of exposure, constrained by the following rates: sublimation of ice, surface erosion of debris, and accumulation of supraglacial debris (Eqs. 14 and 15). If we assume a constant inheritance in the englacial debris, it is only possible to fit the model to the set of samples that show monotonically decreasing nuclide concentrations with increasing depth. The model is therefore fitted to (i) all measured nuclide samples in the supraglacial debris, S1, and (ii) samples of englacial debris, units E1 and E2.

While the recycled paleo-surfaces S2, S3, and S4 are not included in the fitting of the forward model, they are utilized for burial dating to further constrain the age of ice emplacement during optimization of the modeled parameters. The general idea and reasoning for applying burial dating to the recycled paleo-surface samples is as follows. The high cosmogenic nuclide concentrations in the recycled paleo-surfaces are the result of extended periods of exposure of the debris prior to entrainment in the Middle ice. The debris in these samples was part of a surface that was overridden during the latest advance of glacial ice into Ong Valley. Hence, this paleo-surface debris must have been buried at the time the Middle ice was deposited. The burial age obtained from the burial dating of the debris from these recycled paleo-surface samples should then reflect the timing of the Middle ice emplacement. Further, the burial age of these samples cannot display ages that are younger than the event in which they got buried. Therefore, the minimum burial age for any of the recycled paleo-surface sections (S2–4) serves as the maximum age for when the Middle ice was emplaced in Ong Valley.

Burial dating of the recycled paleo-surfaces is applied to their inherited nuclide concentrations. The inherited nuclide concentration is calculated by subtracting the modeled post-emplacement nuclide concentrations from the total measured nuclide concentrations. The apparent burial age of the recycled paleo-surface debris is then determined from the nuclide ratio of the calculated inherited nuclide concentrations. The apparent burial age is the duration of burial inferred from a pair of nuclide measurements under the assumption that a sample has experienced a single period of exposure followed by a single period of burial. In reality, the sample could have experienced multiple shorter periods of burial; however, the calculated apparent burial age is the maximum single period a sample has been buried for. The age of ice emplacement for the Middle ice is then limited by the minimum apparent burial age for any of the recycled paleo-surface samples S2–4. This constraint is incorporated into the model fitting during parameter optimization such that when adding the nuclide concentrations lost by decay during burial, the nuclide ratio does not exceed that of the surface exposure production ratio. This burial constraint is applied to all samples that are not used for forward model fitting, and that make up the recycled surface material S2, S3, and S4.

There are two benefits of applying burial dating to a subset of the samples. First, the burial dating is used to constrain the optimization of the modeled parameters. The concept of burial dating becomes important because of the constraint that, in effect, the samples are not allowed to have a burial age less than zero at the time they are incorporated into the ice. Further, any given sample within the ice core cannot have been buried for less time than the deposit of the ice that encloses it, hence this burial constraint will provide us with a maximum depositional age of the Middle ice.

Secondly, we apply burial dating to determine the burial age of the samples. That is, after we have identified a best-fitting model for the nuclide concentrations produced after ice emplacement, we compute apparent burial ages for samples from the recycled paleo-surface units. This is done to evaluate whether there is a general agreement between the burial age of the samples or if there is variation in the ages which would indicate a more complex history or a variable source of the debris.

6.3 Model results

By fitting the forward model prediction to measured nuclide concentrations from the englacial debris sample (S1, E1, and E2) and applying the burial constraint to sections of recycled paleo-surface debris (S2, S3, and S4), we are able to constrain the age of ice emplacement, sublimation rate, and surface erosion rate for the Ong Valley Middle ice. The results of a 10 000-iteration Monte Carlo simulation provide an ice emplacement age of 2.95 + 0.18 $/$ −0.22 Ma for the Middle ice, with a best-fit χ² of 3.75 + 0.98 $/$ −0.45. The best-fitting sublimation rate since emplacement is 22.86 + 0.10 $/$ −0.09 m Myr⁻¹, with a surface erosion rate of 0.206 + 0.013 $/$ −0.017 m Myr⁻¹. The results of the simulation are not normally distributed, and the best-fit values are therefore reported as the 50th percentile with error bounds given by the 16th and 84th percentile (Figs. 8 and S2). The inherited nuclide concentration for ¹⁰Be, ²¹Ne, and ²⁶Al (the initial nuclide concentration present in the ice mass at the time of deposition 2.95 Myr ago) is 0.101 + 0.018 $/$ −0.017 × 10⁶, 9.2 + 1.4 $/$ −4.4 × 10⁶, and 0.8249 + 0.0062 $/$ −0.0031 × 10⁶ atoms g⁻¹, respectively.

Figure 8 (a–c) Cumulative distribution of 10 000 fitted Monte Carlo simulation results showing the 50th percentile (red line) with error bounds given by the 16th and 84th (pink lines). (f–e) Density plot as paired distribution of Monte Carlo simulation results separated into 1000 bins with yellow being high density and blue low density. Both the vertical and horizontal axes (d–f) are truncated to the 98 % (0.01 and 0.99) confidence interval.

The measured debris concentrations in the ice segments range between 0 for clean ice and 0.57 for debris-rich ice, with an average of 0.13 by weight for the core. The average debris concentration for the best fit in the sublimated ice which produced the supraglacial debris over a period of 2.95 Myr is 0.036–0.034 by weight and therefore an average ice–debris density of ∼ 0.931 g cm⁻³. Further, the resulting sublimation rate of 22.86 m Myr⁻¹ over the span of 2.95 Myr, results in a total ice surface lowering of 67.6 m.

The model predicts nuclide concentrations at depth similar to those measured in the supraglacial debris (S1) and englacial debris (E1–2), and it is therefore evident that these units can be explained with the exposure model (Fig. 9a–c). Further, the paired-nuclide plot (Fig. 9d–f) clearly shows the distinction between the englacial debris data set (S1, E1–2) explainable by the model and the paleo-surface samples (S2–4) having high nuclide concentrations and low paired-nuclide ratio and hence different origins which require the addition of a complex burial and exposure history.

Conceptually, higher nuclide concentrations require long surface exposure, and any disequilibrium in the paired-nuclide ratio (ratio below the steady-state-erosion zone) is the effect of burial after exposure (Lal, 1991). However, the presence of a significant inherited nuclide inventory could result in surface and subsurface samples having ratios below the production ratio, therefore indicating an apparent burial age. In Fig. 9d–f this is observed as both the predicted and all measured nuclide concentration ratios fall within the burial zone and not near the steady-state-erosion zone.

Figure 9 Measured and modeled cosmogenic nuclide concentrations as (a–c) mass depth below the surface for (a) ¹⁰Be, (b) ²¹Ne, and (c) ²⁶Al and (d–f) paired-nuclide exposure–burial diagram for (d) ²⁶Al–¹⁰Be, (e) ¹⁰Be–²¹Ne, and (f) ²⁶Al–²¹Ne pairs. Blue lines indicate the model predicted cosmogenic nuclide concentration from the pit surface to the bottom of the ice core. Colored boxes (a–c) and sample colors (d–f) indicate debris source as detailed in Fig. 7. In (d)–(f), solid black lines show the steady-state erosion zone, and dashed lines show the burial lines as million-year decay isochrons (see Sect. 3 for more details). The measured nuclide concentration for each sample is represented by a shaded ellipse of its 1σ uncertainty. The black line connects the sample ellipses from the surface of the pit down to the bottom of the ice core. ^* signifies nuclide concentrations normalized to respective surface production rate.

6.4 Minimum exposure age

We find the absolute minimum exposure age of the Middle ice to be 1.331 + 0.020 $/$ −0.024 Ma, with a sublimation rate of 24.70 + 0.71 $/$ −0.56 m Myr⁻¹. This age is derived from the minimum possible number of assumptions about the geologic history of the samples. For a surface sample, the apparent age is the calculated age from the measured nuclide concentration assuming a sample has experienced a single event of exposure, zero surface erosion, and no burial during that time period. Under such assumptions, a surface sample's apparent exposure age serves as the minimum exposure age. Therefore, the minimum age for the ice emplacement is obtained using the assumption that the inherited nuclide concentration for ¹⁰Be, ²¹Ne, and ²⁶Al is equal to the minimum concentrations measured throughout the core (0.12 × 10⁶, 6.42 × 10⁶, and 0.41 × 10⁶ atoms g⁻¹, respectively) with zero surface erosion.

6.5 Burial dating of paleo-surface debris

As evident in Fig. 9, the paleo-surface samples have elevated nuclide concentrations and do not fit our modeled predictions. There is no scenario in which these samples can be explained solely by our forward exposure model which includes only sublimation and erosion. Therefore, these samples must have experienced significant periods of surface exposure prior to subglacial entrainment. Further, in order to have a lower paired-nuclide ratio than predicted (Fig. 9d–f), the samples must have experienced at least one period of burial. Hence, these observations were the reasons for the inclusion of burial dating in our model.

Similar to the burial dating constraint added to the forward model (Sect. 6.2), we can determine the burial age of these paleo-surface samples by first subtracting modeled post-depositional nuclide concentrations at the sample depths from the measured concentrations. This yields an estimate of the nuclide concentrations present in the paleo-surface samples (S2–4) at the time they were buried less the effect of subsequent radioactive decay. The choice to only fit the model to a subset of samples is based on the assumption that the paleo-surface samples have different geological history and thus different nuclide inheritance from that of the englacial debris (E1–2) samples. Therefore, the estimated inherited nuclide concentration for these paleo-surfaces obtained from this subtraction is different from the inherited nuclide concentrations inferred from the model fitting. From these inherited nuclide concentrations in the paleo-surface samples, we can then solve for the burial age which would cause a sample exposed at the surface (plotting on the simple exposure line) to have paired-nuclide ratios as shown in Fig. 9. Uncertainties of the burial ages are derived from the same Monte Carlo simulation used to generate uncertainty estimates for the model parameters.

The burial ages for the three paired-nuclide ratios, ²⁶Al $/$ ¹⁰Be, ¹⁰Be $/$ ²¹Ne, and ²⁶Al $/$ ²¹Ne for S2, are 3.21 ± 0.20, 4.20 ± 0.27, and 3.69 ± 0.21 Ma, respectively. The paleo-surface S3 and S4 indicate longer periods of burial, with S3 having burial ages of 4.33 ± 1.00, 7.58 ± 0.61, and 6.24 ± 1.35 Ma and S4 having burial ages of 5.06 ± 0.25, 6.61 ± 0.12, and 5.78 ± 0.15 Ma, respectively, for each of the three nuclide pairs. Figure 10 shows the paired nuclide ratios for each of the paleo-surfaces as their apparent burial ages.

Figure 10 Paired-nuclide diagram for (a) ²⁶Al–¹⁰Be, (b) ¹⁰Be–²¹Ne, and (c) ²⁶Al–²¹Ne pairs. Solid black lines show the steady-state erosion zone, and dashed lines show the burial lines as million-year decay isochrons. The nuclide concentration for each data point is represented by a shaded ellipse of its 1σ uncertainty. Grey data points show the measured nuclide concentrations for the paleo-surface samples as described in Fig. 9d–e. Shaded red data ellipses show the resulting nuclide concentrations when subtracting the modeled nuclide concentration from the measured. The burial age (dashed isochron lines) for which a sample lies represents the apparent burial age and is the maximum single period the sample has been buried since current time. Color shades refer to the different paleo-surface units as described in Fig. 7. ^* signifies nuclide concentrations normalized to respective surface production rate.

7 Discussion

7.1 Sublimation rate

The sublimation rate is tightly constrained between 22.77 and 22.96 m Myr⁻¹ (Fig. 8b) and independent of the age and erosion rate (Fig. 8d, e). With increasing sublimation rate, a sample having low nuclide concentration at deeper depth (caused by increased shielding mass) approaches the ice surface more rapidly. Having spent less time near the ice surface, a sample found in the top meter of the ice will have a much lower total nuclide concentration than a sample at the lower part of the supraglacial debris. This difference in nuclide concentration between the uppermost ice sample and the bottom supraglacial debris sample allows for the sublimation rate to be well constrained.

Previous estimates of the sublimation rate in Ong Valley range between 19 and 22 m Myr⁻¹ (Bibby et al., 2016), whereas sublimation rates of buried ice masses determined using cosmogenic nuclides elsewhere in the Transantarctic Mountains range between 0.7 and 37 m Myr⁻¹ (Schäfer et al., 2000; Ng et al., 2005; Morgan et al., 2010b). These rates broadly agree with a sublimation rate of 22.86 m Myr⁻¹ as reported here in Ong Valley. Sublimation rates obtained from modeled water vapor diffusion are orders of magnitude higher: 100–500 m Myr⁻¹ (Hindmarsh et al., 1998; Mckay et al., 1998; Kowalewski et al., 2006; Hagedorn et al., 2007; Mckay, 2008; Schorghofer, 2009). Such higher sublimation rates by orders of magnitude would suggest a total surface elevation lowering of 300–1500 m as compared to our calculated total ice surface lowering of ∼ 68 m and are inconsistent with glacial moraine elevations and field observations in Ong Valley.

The sublimation rate represents an average rate since the ice emplacement ∼ 3 Myr ago. Most likely the sublimation has been decreasing over time as the supraglacial debris thickens (Mackay and Marchant, 2016). However, this relationship is uncertain, and therefore we do not account for it.

7.2 Erosion rate

The erosion rate and age of the ice are well constrained within an erosion–age tradeoff (Fig. 8f). For an eroding surface, debris of low nuclide concentrations from below will approach the surface at a rate of erosion. With increased surface erosion rate, an older exposure age is required in order to account for the loss of high surficial nuclide concentrations, leading to an expected erosion–age tradeoff.

The majority of Antarctic studies of subaerial surface erosion rate using cosmogenic nuclides are obtained from boulders and bedrock of various lithologies (Marrero et al., 2018, and references therein). Only a few erosion rates have been determined from surficial regolith (Putkonen et al., 2008; Morgan et al., 2010a; Bibby et al., 2016). While Bibby et al. (2016) found a 0.89 m Myr⁻¹ for the Middle drift, a range between 0.2 and 2 m Myr⁻¹ has been observed in McMurdo Dry Valleys (Putkonen et al., 2008; Morgan et al., 2010a). Therefore, an erosion rate of 0.206 m Myr⁻¹ as reported here for the supraglacial debris is in agreement.

7.3 Englacial debris concentration

In Ong Valley, we measured an average debris concentration of 0.13 by weight in the ice (Eq. 1) which is in the same range as measurements made in Beacon Valley (0.085 by Marchant et al., 2002; Morgan et al., 2010a). While the modeled debris concentration of 0.035 by weight in the sublimated ice over a span of 2.95 Myr is lower than measured debris concentration of buried ice in Antarctica, it is consistent with the expectation that the debris content increases towards the bottom of glacial ice due to subglacial entrainment of the debris. Thus, it is expected that the modeled debris concentration for the sublimated ice here results in a lower concentration than measured in the remaining basal ice.

7.4 Mixing layer

Predicted cosmogenic nuclide concentrations in the supraglacial debris decrease with depth at a higher rate than measured nuclide concentrations (Fig. 9a–c). This leads to a systematic misfit between observations and model predictions. By either decreasing the sublimation rate, increasing the erosion rate, and/or decreasing the age of ice emplacement, a steeper predicted cosmogenic nuclide depth profile can be obtained for the supraglacial debris. However, neither of these scenarios will result in a better fit for the near-surface pit samples. The difficulty of fitting the forward model to the near-surface pit samples suggests that partial vertical mixing of the supraglacial debris may have occurred.

As the supraglacial debris is accumulating due to sublimation, debris having low cosmogenic nuclide concentrations from below will mix with debris of higher surficial cosmogenic nuclide concentrations. Therefore, any (partial or full) vertical mixing of the supraglacial debris would cause a decrease in the cosmogenic nuclide inventory in any above-lying sample. Without accounting for any vertical mixing of the supraglacial debris, the model predictions result in an overestimation of the cosmogenic nuclide concentration for the near-surface pit samples and an underestimation of newly accumulating supraglacial debris from the ice surface as observed in Fig. 9a–c. Further, any vertical mixing would decrease the nuclide ratio, which would explain why all paired-nuclide ratios for the supraglacial debris samples plot below the steady-state erosion zone (Fig. 9d–f).

It has previously been suggested that no vertical mixing occurs in the supraglacial debris layers in Ong valley (Bibby et al., 2016) and supraglacial debris layers studied in Beacon Valley (Morgan et al., 2010a, b). While the current measured nuclide profile may not reflect a fully mixed zone as seen elsewhere in temperate climate with bioturbation (Perg et al., 2002), a partially mixed supraglacial debris layer is likely the result of active polygon formation found at the surface of the Middle drift.

7.5 Exposure ages from boulders

In general, a boulder having experienced a single period of exposure that is equal to the ice emplacement age of the Middle drift should display concordant ¹⁰Be and ²¹Ne ages. However, all measured boulders show apparent ¹⁰Be exposure ages younger than those of ²¹Ne (Table 2) and are therefore inconsistent with a simple exposure having negligible erosion.

From the ¹⁰Be $/$ ²¹Ne ratio (Fig. 11) it is observed that the surface boulders have experienced erosion while exposed to cosmic rays at the surface as the paired-nuclide ratio lays within the steady-state erosion zone (see details in Sect. 3). Three outliers of the east lateral moraine boulders (214, 217, and 218; Fig. 11a) show neither age nor erosion rates that agree with continuous exposure and lie below the steady erosion zone in a region of intermittent exposure. Thus, these boulders show a complex exposure history having experienced at least one period of burial at some point in time.

Figure 11 ¹⁰Be–²¹Ne paired-nuclide diagram of the boulder samples. Solid black lines are the simple exposure line and steady-state erosion line, which marks the zone of continuous exposure. Blue lines are isolines of constant steady erosion (cm Myr⁻¹), and dashed black lines are isolines of constant exposure age (Myr). The measured nuclide concentration for each sample is represented by a red dot with red shading of its 1σ uncertainty. ^* signifies nuclide concentrations normalized to respective surface production rate.

A more realistic exposure age and erosion rate can be determined for boulders having a nuclide ratio within the steady-state erosion zone. By assuming a single period of continuous exposure at a steady-state erosion, we can solve for both the exposure and erosion rates as detailed in Balco et al. (2014). The results of a 10 000-iteration Monte Carlo simulation using this procedure are shown in Fig. 12. Some samples permit infinite ages at a steady erosion rate if the ¹⁰Be $/$ ²¹Ne nuclide ratio lies outside of the continuous exposure zone. Therefore, only samples permitting finite age–erosion rate solutions are shown in Fig. 12.

Figure 12 Exposure ages and erosion rates for the Middle drift surface and moraine boulders.

The lateral moraine boulders having a finite age–erosion rate solution display a 68 % confidence bound on the age of 1.0–3.9 Ma and erosion rates of 0.20–0.48 m Myr⁻¹ (Fig. 12a). From field observations, the lateral moraine from which boulder measurements were sampled appears to be a younger recessional moraine and therefore not an indication of a maximum extent for the Middle drift, which is observed at higher elevation. These observations would suggest that boulders from this lateral moraine have most likely been disturbed after ice emplacement. Thus, the ¹⁰Be $/$ ²¹Ne ratio age for the moraine is more likely to represent a minimum limiting age of ice emplacement (Hallet and Putkonen, 1994; Putkonen and Swanson, 2003).

The Middle drift surface boulders have a 68 % confidence bound on the age of 1.8–3.5 Ma with erosion rates ranging between 0.05 and 0.21 m Myr⁻¹ and therefore agree with our modeled age of 2.96 Ma for the Middle drift ice (Fig. 12b). One outlier (11-OV-ER-117) has steady-state erosion for a period greater than several half-lives of ¹⁰Be (≳ 5 Myr) and therefore contains no age information. We attribute this increased nuclide concentration to an extended period of exposure prior to deposition in Ong Valley.

7.6 Multiple glacial events

Samples having a burial age equal to that of the ice emplacement age are considered to have been derived from a paleo-surface exposed in Ong Valley during ice advance of the Middle ice. After entrainment into the advancing ice, the paleo-surface material was immediately buried under a thick layer of ice and shielded from the cosmic-ray flux.

As mentioned in Sect. 6.2, the age of the ice determined from modeled nuclide predictions is constrained such that the ice cannot be older than the minimum burial age obtained from any sample across the three paired-nuclide ratios. We find that the minimum burial age of 3.21 ± 0.20 Ma for S2 agrees with the age of the Middle ice (2.95 + 0.18 $/$ −0.21 Ma). Hence, S2 is likely to have been at the surface during the glacial advance leading to the deposition of the Middle ice at ∼ 3 Ma. S2 is found at depths in between E1 and E2 (Fig. 7), which have no prior exposure history. This would suggest that S2 is not in stratigraphic order and has been mixed up into the ice during advance of the Middle ice into Ong Valley.

Both S3 and S4 display older burial ages which are not uniform across the three paired-nuclide ratios (Fig. 10). This suggests that both S3 and S4 have experienced a complex exposure–burial history of prior entrainment which goes beyond the exposure history of the Middle drift ice. The estimated burial ages represent the minimum total burial time a sample has experienced but also the maximum burial time of a single burial event. For S3 and S4, additional burial time is needed beyond the age of the ice (> 2.95 Ma) and must have experienced multiple periods of burial. The simplest explanation is to assume that during advance of the Middle ice, Ong Valley looked similar to today such that the Middle ice advanced over an already existing ice-cored drift unit. Perhaps this is now preserved as the Old drift (Fig. 2). If this is the case, then S3 and S4 units presumably are debris that was buried for some time in the older ice and then buried again in or under the Middle ice. We use “in or under” because it is possible that (i) we cored through the Middle ice into a stratigraphically underlying mass of older ice. However, it is also possible that (ii) sections of the older ice were entrained during advancement of the Middle ice, and we then have a mixture of older and younger ice. The core did not display an obvious stratigraphic boundary. In fact, most of the observed grain in the ice–debris mass is at steep angles and disturbed (Sect. 4.2), which would tend to favor a mixing of the ice hypothesis. Regardless, either is possible, and geochemical analysis of the ice could potentially help resolve this.

The uncertainty associated with the age of this older, separate underlying ice mass is greater compared to the Middle ice due to the complexity associated with the exposure–burial dating. That is, the burial age obtained here is the apparent burial of a single event. However, a sample could have experienced multiple shorter periods of exposure–burial events which are not accounted for. We therefore only report an estimate with minimum and maximum constraints. A sample that experiences a single period of exposure at the surface has a nuclide concentration ratio dependent on the duration of exposure. When buried to a depth where nuclide production is significantly reduced, the change in ratio is primarily dependent on the radioactive decay associated with the duration of burial. Therefore, when solving for the burial age for each of these paleo-surface units, we can also solve for the exposure age that has occurred prior to burial by the Middle ice (Balco and Rovey, 2008).

With S2 representing a surface sample from a supraglacial debris prior to deposition of the Middle ice, we find the surface exposure ages of S2 to be 0.163 ± 0.058, 0.268 ± 0.046, and 0.268 ± 0.046 Ma for the paired-nuclide ratios ²⁶Al $/$ ¹⁰Be, ¹⁰Be $/$ ²¹Ne, and ²⁶Al $/$ ²¹Ne, respectively. This suggests that the age of an underlying ice mass is at least that of the minimum burial age plus the exposure age for S2. On the contrary, the age of the ice in which a sample is embedded cannot be older than the burial age of a given sample. Then, the maximum burial age across the three paleo-surfaces S2–4 for any of the paired-nuclide ratios must serve as the upper bound for the age. Therefore, there is no scenario in which this separate, older underlying ice mass can be younger than 3.3 Ma or older than 7.6 Ma. However, a more likely age for this older deposit would be 4.3–5.1 Ma, which is the burial age obtained from the ²⁶Al $/$ ¹⁰Be ratio for S3 and S4 as it is unclear whether or not S3 and S4 are the same or different units.

It is difficult to determine whether or not there is any defined boundary between the older and younger ice masses. With an increase in nuclide concentrations downcore and, in addition, samples from E2, S3, and S4 appearing to form mixing arrays in the paired-nuclide diagrams shown in Fig. 9, it appears likely that S3 is a mixture of a high-nuclide-concentration end-member, which may be represented by S4, and a low-nuclide-concentration end-member represented by E2. However, a boundary or transition most likely exists between E2 and S4.

7.7 Antarctica during the Pliocene

The ages reported here coincide with the Pliocene epoch (5.3–2.6 Ma). Research on Pliocene climate and how it affected the Antarctic ice sheet has gained much attention as a likely analog for modern anthropogenic warming (Dolan et al., 2018). During the Pliocene epoch there is evidence of prominent glacial deposits, in which two are identified in the Southern Hemisphere as globally recognizable glaciations (summarized in De Schepper et al., 2014): one occurring during the early Pliocene (ca. 4.9–4.8 Ma) and another during the late Pliocene (ca. 3.3 Ma), also identified as the Marine Isotope Stage (MIS) M2 glaciation. The latter is followed by a warmer than present mid-Piacenzian warming period (mPWP; ∼ 3.3–3.0 Ma) (Haywood et al., 2013; De Schepper et al., 2014; Dowsett et al., 2016). This warming period ends by a late Pliocene cooling, after ∼ 3 Ma, leading to a global glaciation around the Pliocene–Pleistocene boundary (De Schepper et al., 2014).

Because the uncertainty of the ice emplacement age (± 0.2 Ma) exceeds both the 40 and 100 kyr climate cycles of the Pliocene epoch (Caballero-Gill et al., 2019), we are not able to accurately relate the deposition of the Middle ice to an individual glacial event. Furthermore, the age of 2.95 + 0.18 $/$ −0.22 Ma for the Middle ice emplacement, which requires an East Antarctic Ice Sheet (EAIS) elevation greater than 200 m above present, cannot confidently be assigned to either the warmer period prior to 3 Ma or the cooler period after 3 Ma. Balter-Kennedy et al. (2020) concluded that glacial deposits recording a higher than present EAIS elevation at Roberts Massif, a nearby location in the Transantarctic Mountains, most likely postdated the mPWP. Therefore, if the ice advance in Ong Valley was correlative with that at Roberts Massif, it would also be associated with the 3 Ma cooling. However, this is speculative.

The oldest englacial debris that we have dated in Ong Valley is dated at ∼ 4.3–5.1 Myr. Although the dated age range is rather wide due to complexities resulting from old age and exposure–burial dating, it is still direct evidence of an EAIS expansion and local ice expansion during that time. This dated age suggests that the ice sheet expansion predated the MIS M2 cooling event and possibly coincided with the early Pliocene global glaciation (ca. 4.9–4.8 Ma). If in fact the older ice is still present below the Middle ice mass, then it did not melt during a period of warming. Thus, additional evidence indicating whether or not two ice units are present would be important in understanding the climate during the Pliocene epoch. Since ∼ 3 Ma there has never been a comparable ice sheet expansion in Ong Valley as was seen in the early/middle and late Pliocene. The only notable but small ice sheet advance or stagnation evident in Ong Valley is the Young drift dated at 11–13 ka.

8 Conclusions

Glacial ice is a well-known paleoclimate archive. Great efforts have been made to find ice older than 1 Ma since the paucity of ice beyond million years of age creates uncertainties for future climate predictions. In Ong Valley, Antarctica, the Middle drift harbors a large ice mass buried 62 cm below the surface of supraglacial debris. We collected a 944 cm long ice core and measured concentrations of the cosmic-ray-produced nuclides ¹⁰Be, ²⁶Al, and ²¹Ne from the englacial debris and samples from the supraglacial debris directly above it. We developed a numerical forward model which predicts the accumulation of cosmogenic nuclides in the englacial debris and the above-lying supraglacial debris during a single event of exposure, constrained by sublimation, surface erosion, and accumulation of supraglacial debris. The modeled nuclide concentrations are then fitted to the measured nuclide concentrations in the ice core.

A downcore increase in measured nuclide concentrations suggests that sections of englacial debris consist of both subglacially entrained debris and recycled paleo-surfaces having a complex exposure–burial history prior to entrainment. This allows us to apply a combination of exposure and burial dating to the forward model. We find the age of the Middle drift ice mass to be 2.95 Ma, with a constant ice sublimation rate of 22.86 m Myr⁻¹ and surface erosion rate of 0.206 m Myr⁻¹. Cosmogenic nuclide exposure dating of surface boulders belonging to the surface of the coring site are consistent with the modeled age of ∼ 3 Ma for the ice emplacement.

Exposure–burial dating on the englacial paleo-surface debris reveals that the lower section of the ice core belongs to a separate and older deposit, emplaced ∼ 4.3– 5.1 Myr ago. We interpret this lower section as a portion of an older ice mass either in situ or transported during emplacement of the younger ice. The ages of the two separate ice masses found below the Middle drift can be directly related to glacial advances. These findings provide direct evidence of an Antarctic ice sheet that was larger than present during the early and late Pliocene epoch.

Furthermore, we show that exposure–burial dating of cosmogenic nuclides measured in situ in basal ice debris layers can be used for age constraint of past ice advance. Specifically, we have debris layers in one ice core that suggest three different burial ages, in which at least two of them are dated to be older than the age of the ice itself. This is important for understanding in situ cosmogenic nuclide data from out-of-context subglacial sediment.

Collectively our results show that the continental ice sheet advanced into Ong Valley repeatedly, and evidence of at least two such advances at 2.95 and 4.3–5.1 Ma is still preserved in lateral moraines, drifts, and stacked ice masses. Since 2.95 Ma the only evidence of ice advance or stagnation in the Ong Valley was ∼ 10 kyr ago.